//  Copyright (c) 2012 M.A. (Thijs) van den Berg, http://sitmo.com/
//
//  Use, modification and distribution are subject to the MIT Software License. 
// (See accompanying file LICENSE.txt or copy at http://www.stdfin.org/LICENSE.txt)
// 

#ifndef STDFIN_DISTRIBUTION_NORMAL_HPP
#define STDFIN_DISTRIBUTION_NORMAL_HPP
#include <cmath>

#ifdef USE_BOOST
    #include <boost/math/distributions/normal.hpp>
#endif

namespace stdfin {

#ifdef USE_BOOST

    // -------------------------------------------------------------
    // USE BOOST
    // -------------------------------------------------------------
    inline double normal_cdf(double x)
    {
        return boost::math::cdf( boost::math::normal(), x);
    }
    
    inline double normal_pdf(double x)
    {
        return boost::math::pdf( boost::math::normal(), x);
    }

#else
    // -------------------------------------------------------------
    // DON'T USE BOOST
    // -------------------------------------------------------------
    
    // This function is accurate to appro 10 decimal places
    // Based on the paper "BETTER APPROXIMATIONS TO CUMULATIVE NORMAL FUNCTIONS" Graeme West
    // orginal source: Hart, J. (1968), Computer Approximations, Wiley. Algorithm 5666 for the error function.
    inline double normal_cdf(double x)
    {
        double Cumnorm;
        double build;
    
        double XAbs = std::abs(x);
        if (XAbs > 37.0) 
            Cumnorm = 0.0;
        else {  
            double Exponential = exp( -XAbs*XAbs / 2.0);
            if (XAbs < 7.07106781186547) { 
                build = 3.52624965998911E-02 * XAbs + 0.700383064443688;
                build = build * XAbs + 6.37396220353165;
                build = build * XAbs + 33.912866078383;
                build = build * XAbs + 112.079291497871;
                build = build * XAbs + 221.213596169931;
                build = build * XAbs + 220.206867912376;
                Cumnorm = Exponential * build;
                build = 8.83883476483184E-02 * XAbs + 1.75566716318264;
                build = build * XAbs + 16.064177579207;
                build = build * XAbs + 86.7807322029461;
                build = build * XAbs + 296.564248779674;
                build = build * XAbs + 637.333633378831;
                build = build * XAbs + 793.826512519948;
                build = build * XAbs + 440.413735824752;
                Cumnorm = Cumnorm / build;
            } else {
                build = XAbs + 0.65;
                build = XAbs + 4 / build;
                build = XAbs + 3 / build;
                build = XAbs + 2 / build;
                build = XAbs + 1 / build;
                Cumnorm = Exponential / build / 2.506628274631;
            }
        }
        if (x>0) {
            Cumnorm = 1.0 - Cumnorm;
        }
        return Cumnorm;
    }
    
    inline double normal_pdf(double x)
    {
        return exp(-0.5*x*x) * 0.3989422804014326779399460599343818684758586311649346; 
    }
    
#endif


} //  namespace

#endif // include guard